In [1]:
import sympy
import sympy.physics.mechanics as mech
import scipy.integrate
import numpy as np
import matplotlib.pyplot as plt
import collections
import os
import sympy_utils
from sympy.printing.theanocode import theano_function
sympy.init_printing()
sympy.physics.vector.init_vprinting()
%matplotlib inline
%load_ext autoreload
%autoreload 2
plt.rcParams['lines.markersize']= 10
plt.rcParams['font.size']= 12
In [2]:
phi, theta, psi, P, Q, R, V, alpha, beta = \
sympy.symbols('phi, theta, psi, P, Q, R, V, alpha, beta')
eta, eta_dot, alpha_l, alpha, beta_l, beta, s, t = \
sympy.symbols('eta, eta_dot, alpha_l, alpha, beta_l, '
'beta, s, t')
states_sub = {phi(t): phi, theta(t): theta,
P(t): P, Q(t): Q, R(t): R,
V(t): V, alpha(t): alpha, beta(t): beta,
s(t): s,
eta(t).diff(t): eta_dot, eta(t): eta}
In [3]:
frame_i = mech.ReferenceFrame('i')
frame_b = frame_i.orientnew('b', 'Body', (psi(t), theta(t), phi(t)), '321')
frame_b.set_ang_vel(frame_i, P(t)*frame_b.x + Q(t)*frame_b.y + R(t)*frame_b.z)
frame_w = frame_b.orientnew('w', 'Body', (beta(t), -alpha(t), 0), '321')
frame_f = frame_b.orientnew('f', 'Axis', (eta(t), frame_b.x))
frame_l = frame_f.orientnew('l', 'Body', (-beta_l(t), -alpha_l(t), 0), '321')
In [4]:
point_cm = mech.Point('cm')
point_cm.set_vel(frame_i, V*frame_w.x)
point_cm.set_vel(frame_b, 0)
point_p = point_cm.locatenew('p', -frame_f.y*s)
point_p.set_vel(frame_f, 0)
In [5]:
vel_p = point_p.v1pt_theory(point_cm, frame_i, frame_f)
vel_p
Out[5]:
In [6]:
alpha_l_expr = sympy.atan(vel_p.dot(frame_f.z)/ vel_p.dot(frame_f.x)).subs(states_sub)
alpha_l_expr
Out[6]:
In [7]:
beta_l_expr = sympy.atan(vel_p.dot(frame_f.y)/ vel_p.dot(frame_f.x)).subs(states_sub)
beta_l_expr
Out[7]:
In [8]:
C_L_alpha = sympy.Symbol('C_L_alpha')
C_D_0 = sympy.Symbol('C_D_0')
C_L_DM = sympy.Symbol('C_L_DM')
C_L_0 = sympy.Symbol('C_L_0')
c_bar = sympy.Symbol('c_bar')
In [9]:
C_L = C_L_alpha*alpha_l(t) + C_L_0
C_L
Out[9]:
In [10]:
k_CD_CL, rho, S = sympy.symbols('k_CD_CL, rho, S')
C_D = k_CD_CL*(C_L - C_L_DM)**2 + C_D_0
C_D
Out[10]:
In [11]:
q = rho*vel_p.dot(vel_p)/2
q
Out[11]:
In [12]:
L = -C_L*q*c_bar*frame_l.z
L
Out[12]:
In [13]:
D = -C_D*q*c_bar*frame_l.x
D
Out[13]:
In [14]:
F_A = (L + D).subs(states_sub)
F_A_b = F_A.to_matrix(frame_b).subs({alpha_l(t): alpha_l_expr, beta_l(t): beta_l_expr})
F_A
Out[14]:
In [15]:
x_ac = sympy.Symbol('x_ac')
M_A = (-x_ac*frame_b.x + s*frame_f.y).cross(F_A).subs(states_sub)
M_A_b = M_A.to_matrix(frame_b).subs({alpha_l(t): alpha_l_expr, beta_l(t): beta_l_expr})
M_A
Out[15]:
In [16]:
def gen_aero_force_moment_funcs(F_A_b, M_A_b, const):
# generate a fast function to find the aero force
f_aero = theano_function(
sympy.symbols('P Q R V alpha beta eta eta_dot s'),
[sympy.Matrix([F_A_b, M_A_b]).subs(states_sub).subs(const)],
on_unused_input='warn', allow_input_downcast=True)
return f_aero
In [17]:
def compute_aero(f_aero, P, Q, R, V, alpha, beta, eta_l, eta_dot_l, eta_r, eta_dot_r, const):
"""
Given state of aircraft, computes instantaenous aerodynamic force.
"""
if abs(V) < 1e-3:
return np.zeros(3), np.zeros(3)
span = const['span']
ode_left = scipy.integrate.ode(lambda s: f_aero(P, Q, R, V, alpha, beta, eta_l, eta_dot_l, s))
#ode_left.set_integrator('dopri5')
ode_left.set_initial_value([0,0,0,0,0,0], 0)
ode_left.integrate(span/2.0)
FM_left = ode_left.y
ode_right = scipy.integrate.ode(lambda s: f_aero(P, Q, R, V, alpha, beta, -eta_r, -eta_dot_r, -s))
#ode_right.set_integrator('dopri5')
ode_right.set_initial_value([0,0,0,0,0,0], 0)
ode_right.integrate(span/2.0)
FM_right = ode_right.y
FM = FM_left + FM_right
F = FM[:3]
M = FM[3:]
return F, M
In [18]:
def calc_aero_forces(f_aero, t, eta_l, eta_r, P, Q, R, V, alpha, beta, const):
"""
Computes force history due to flapping motion of wing.
"""
# do simulation
n = len(t)
F_b_vect = np.zeros((n,3))
M_b_vect = np.zeros((n,3))
dt = np.diff(t)
eta_dot_l = np.diff(eta_l)/dt
eta_dot_l = np.hstack([eta_dot_l[0], eta_dot_l])
eta_dot_r = np.diff(eta_r)/dt
eta_dot_r = np.hstack([eta_dot_r[0], eta_dot_r])
for i in range(n):
F_b, M_b = compute_aero(
f_aero=f_aero, P=P, Q=Q, R=R,
V=V, alpha=alpha, beta=beta,
eta_l=eta_l[i], eta_dot_l=eta_dot_l[i],
eta_r=eta_r[i], eta_dot_r=eta_dot_r[i],
const=const)
F_b_vect[i,:] = F_b
M_b_vect[i,:] = M_b
return F_b_vect, M_b_vect
In [19]:
def plot_aero_force_moment(t, eta_l, eta_r, F_b, M_b):
"""
Plots aerodynamic force history.
"""
plt.figure(figsize=(10,10))
plt.subplot(311)
plt.plot(t, eta_l,'b-')
plt.plot(t, eta_r,'r--')
#plt.xlabel('t, sec')
plt.legend(['left','right'], loc='best')
plt.ylabel(r'$\eta$')
plt.title('wing movement')
plt.grid()
lift = -F_b[:,2]
side = -F_b[:,1]
thrust = F_b[:,0]
plt.subplot(312)
newton_to_grams = 1.0e3/9.8
plt.plot(t, lift*newton_to_grams)
plt.plot(t, thrust*newton_to_grams)
plt.plot(t, side*newton_to_grams)
plt.legend(['lift', 'thrust', 'side'], loc='best')
#plt.xlabel('t, sec')
plt.ylabel('gm')
plt.title('aerodynamic forces')
plt.grid()
plt.subplot(313)
newton_m_to_gram_cm = 1.0e4/9.8
plt.plot(t, M_b*newton_m_to_gram_cm)
plt.legend(['roll', 'pitch', 'yaw'], loc='best')
plt.xlabel('t, sec')
plt.ylabel('gm-cm')
plt.title('aerodynamic moments')
plt.grid()
In [20]:
def show_mean_force_moment(F_b, M_b, alpha):
newtons_to_gm = 1000/9.8
newton_m_to_gm_cm = 1.0e4/9.8
F_b_mean = F_b.mean(0)*newtons_to_gm
M_b_avg_gm_cm = M_b.mean(0)*newton_m_to_gm_cm
lift_avg_gm = -F_b_mean[2]
side_avg_gm = F_b_mean[1]
thrust_avg_gm = F_b_mean[0]
print 'lift:\t{:10.2f} gm\nside:\t{:10.2f} gm\nthrust:\t{:10.2f} gm'.format(
lift_avg_gm, side_avg_gm, thrust_avg_gm)
print 'roll:\t{:10.2f} gm-cm\npitch:\t{:10.2f} gm-cm\nyaw:\t{:10.2f} gm-cm'.format(
M_b_avg_gm_cm[0], M_b_avg_gm_cm[1], M_b_avg_gm_cm[2])
F_fwd = -lift_avg_gm*np.sin(alpha) + thrust_avg_gm*np.cos(alpha)
F_up = lift_avg_gm*np.cos(alpha) + thrust_avg_gm*np.sin(alpha)
print 'fwd: {:f} gm, up: {:f} gm'.format(F_fwd, F_up)
In [21]:
s1_aero_fit = sympy_utils.load_repr(os.path.join('save', 's1_aero_fit.repr'))
const = dict(s1_aero_fit['glide'])
const.update({'span': 0.6, 'x_ac': 0.1,
'c_bar': 0.1, 'S': 1, 'rho': 1.225})
# generate a fast function to find the aero force/moment
f_aero = gen_aero_force_moment_funcs(F_A_b, M_A_b, const)
pi = np.pi
In [22]:
def sim_flap(V, alpha, freq, center, delta, flap_mag, f_aero):
pi = np.pi
beta = 0
t = np.linspace(0, 1.0/freq, 100)
eta_l = flap_mag*np.sin(2*pi*freq*t) + center + delta
eta_r = flap_mag*np.sin(2*pi*freq*t) + center - delta
F_b, M_b = calc_aero_forces(
f_aero=f_aero,
t=t, eta_l=eta_l, eta_r=eta_r,
P=0, Q=0, R=0, V=V, alpha=alpha, beta=beta, const=const)
plot_aero_force_moment(t, eta_l, eta_r, F_b, M_b)
show_mean_force_moment(F_b, M_b, alpha)
In [23]:
sim_flap(V=1, alpha=20*pi/180, freq=5, center=6.5*pi/180, delta=0*pi/180, flap_mag = 100*pi/180, f_aero=f_aero)
In [24]:
sim_flap(V=1, alpha=20*pi/180, freq=5, center=50*pi/180, delta=0*pi/180, flap_mag = 100*pi/180, f_aero=f_aero)
In [25]:
sim_flap(V=1, alpha=20*pi/180, freq=5, center=-50*pi/180, delta=0*pi/180, flap_mag = 100*pi/180, f_aero=f_aero)
In [26]:
sim_flap(V=1, alpha=20*pi/180, freq=5, center=6.5*pi/180, delta=-30*pi/180, flap_mag = 100*pi/180, f_aero=f_aero)
In [27]:
sim_flap(V=1, alpha=20*pi/180, freq=5, center=6.5*pi/180, delta=30*pi/180, flap_mag = 100*pi/180, f_aero=f_aero)
In [28]:
f_aero(*np.ones(9).astype(float))
Out[28]:
In [29]:
f_aero(*np.ones(9).astype(float))
Out[29]:
In [30]:
%%timeit
f_aero(*np.ones(9).astype(float))
In [31]:
compute_aero(f_aero=f_aero,
P=0, Q=0, R=0,
alpha=0.1, beta=0, V=1,
eta_l=0, eta_dot_l=0,
eta_r=0, eta_dot_r=0, const=const)
Out[31]:
In [32]:
%%timeit
compute_aero(f_aero=f_aero,
P=0, Q=0, R=0,
alpha=0.1, beta=0, V=1,
eta_l=0, eta_dot_l=0,
eta_r=0, eta_dot_r=0, const=const)